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We demonstrate why for a sheared gas of hard spheres, described by the SLLOD equations with 
an iso-kinetic Gaussian thermostat in between collisions, deviations of the conjugate pairing rule for 
the Lyapunov spectrum are to be expected, employing a previous result that for a large number of 
particles N, the iso-kinetic Gaussian thermostat is equivalent to a constant friction thermostat, up 
to 1/yN fluctuations. We also show that these deviations are at most of the order of the fourth 
power in the shear rate. 



PACS Numbers: 05.20.-y, 05.45.-a, 05.60.Cd 



The SLLOD equations of motion, combined with Lees- 
Edwards boundary condition Q, were originally pro- 
posed in Refs. and since then they have been con- 
venient tools to calculate the shear viscosity of gases in 
the bulk by means of non-equilibrium molecular dynam- 
ics simulations for many years. These studies consider 
systems with a large number of mutually interacting par- 
ticles that are driven by an external shear rate 7 ■ In 
these studies, the iso-kinetic Gaussian thermostat is an 
artificial way to continuously remove the energy gener- 
ated inside the system due to the work done on it by the 
external shear field, such that a non-equilibrium steady 
state, homogeneous in space, can be reached. The Lya- 
punov spectrum of such systems is of interest since it has 
been shown that the shear viscosity can be related to the 
spectrum [f|[|, which can be numerically obtained as a 
function of the shear rate @ . The analysis of the simula- 
tion data H indicated that the sum of the largest and the 
smallest, the sum of the second largest and the second 
smallest and so on, were the same. The phenomenon of 
such pairing of the Lyapunov exponents is known as the 
Conjugate Pairing Rule, or the CPR. Based on this obser- 
vation, an attempt to prove an exact CPR was made for 
arbitrary inter-particle potentials and arbitrary 7 
and later studies and better simulation techniques Q '. 
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indicated that for systems obeying the SLLOD equations 
of motion, the CPR is not satisfied exactly under these 
general conditions |T^ |. However, any conclusive theoret- 
ical proof regarding the status of an approximate CPR 
for systems under SLLOD equations of motion is absent 
in the literature till now, leaving the problem open for a 
long time. 

The SLLOD equations of motion describe the dynam- 
ics of a collection of N particles constituting a fluid with 
a macroscopic velocity field u(r) = 72/x. For particles of 
unit mass, the equations of motion of the z-th particle, 
in terms of its position and peculiar momentum pj, is 



given by 



pi = Fi - jpiyi - apt 



(1) 



where is the force on the z-th particle due to the other 
particles in the system. The value of a, the coefficient of 
friction representing the iso-kinetic Gaussian thermostat, 
is chosen such that the total peculiar kinetic energy of the 
system, X^P?/^! is a constant of motion in between col- 
lisions. In terms of the positions and the laboratory 
momenta of the particles, Eq. (Q) reads 

ii = v; , Vj = Ft + a-fyiic - av; . (2) 

In the present context, the gas particles are hard 
spheres, which for simplicity are assumed to have unit 
radius. The dynamics of the gas particles consists of 
an alternating sequence of flight segments and instanta- 
neous binary collisions. During a flight, the dynamics of 
the gas particles is described by Eqs. (j^) with Fj = 0. At 
an instantaneous collision between the i-th and the j-th 
sphere, the post-collisional positions and laboratory mo- 
menta (+ subscripts) are related to their pre-collisional 
values (— subscripts) by 



- {(Vj_ - Vj_) • Ay} fiy and 
+ {(Vj_ - Vj-_) • fiy , 



(3) 



while the positions and the velocities of the rest of the 
spheres remain unchanged. Here, fiy is the unit vector 
along the line joining the center of the z-th sphere to the 
j-th sphere at the instant of collision. Note that because 
we applied the iso-kinetic Gaussian thermostat only be- 
tween collisions |l3| , the peculiar kinetic energy changes 
in individual collisions. These changes are random, both 
in magnitude and sign, due to the randomness of the col- 
lision parameters, and hence it is quite likely that the 
system would reach a steady state, where the average 
change of peculiar kinetic energy would be zero. 
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In terms of the 3iV-dimensional vectors R 



(ri,r 2 , ...,r N ), V = 
entry is given by N: 
Eqs. 



3) can be compacted to 
R = V , V = cryCR 
during a flight segment and 



(vi, v 2 , vjy) and N^, whose l-th 
= {5 lti -8 l>j )i Hj /j2{l = l,2,..,N), 



aV 



(4) 



R4 



R 



V+ = V_ - 2 (V_ • N y ) Nj 



at a collision between the «-th and the j-th sphere p4[ . 
Here, C is a 3^x3^ matrix with NxN entries, each of 
which is a 3x3 matrix. In terms of the entry index (I, m), 
in the xyz-basis, Q m = c <5; m (/, m = 1,2,.., N) and 
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Having described the dynamics of the infinitesimal de- 
viation £X = (SR, 5~V) between two typical trajectories 
in the 6-/V-dimensional phase space for a time t as 



SX(t) = L(t) <5X(0) , 



(5) 



the Lyapunov exponents for this system are the loga- 
rithms of the eigenvalues of the matrix A, defined by 



A = lim 

t — >oo 



L(t) 



1/(20 



where L(i) = [L(i)] T L(i). 

It can be shown ]l5| ] that the sufficient condition for 
the CPR to hold exactly for a dynamical system obeying 
Eq. (J5J) is the existence of a constant non-singular matrix 
K satisfying K 2 cxl, such that 



[L(i)] T K L(t) = fiK. 



(6) 



Here, /1 is a scalar function of t. If L(t) satisfies Eq. (g), 
then we call L(t) to be "generalized /i-symplectic" . It is 
easy to show from Eq. (g|) that if L is an eigenvalue of 
L(t), then so is /i 2 /Z; from which the (exact) CPR fol- 
lows. For the situations where the CPR has been proved 
to be exact |]l^ , ^6| -[l8| , only the ^-symplecticity case of 
Eq. (0) (i.e., K = J, where J is the usual symplectic ma- 
trix) has been exploited. In this context, we note that 
despite the similarity between the present problem and 
the one discussed in Ref. H], the elaborate formalism 
developed therein is not applicable here. 

A significant simplification can be achieved by noticing 
that the coefficient of friction a, in the non-equilibrium 
steady state, fluctuates with l/>/~N fluctuations around a 
fixed value ao in the thermodynamic limit p9[ . Thus, to 
calculate the Lyapunov exponents for large N, to which 
we confine ourselves henceforth, a can be replaced by ao 
in Eq. (0), except for a beginning transient time. On 



average, for small 7, a oc j 2 and so is ao- Higher order 
corrections play a role for larger shear rates. 

In the following analysis, we first explore the status of 
the CPR when the coefficient of friction is a constant, ao, 
and then return to the case where the coefficient of fric- 
tion represents an iso- kinetic Gaussian thermostat. The 
detailed derivation of the following results is given else- 
where Jl5| . At present, we focus only on the main points. 

Once ao replaces a in Eq. (Q), we find that in the 
time evolution of (5X over a collision-less flight segment 
between t and t + r is given by 



<5X(t + r) = H(r) £X(t) . 



(7) 



H(r) can be decomposed into 3iV x 3A^ sub- matrices as 



H(r) 



hW(r) hl 2 ](r) 
hl 3 l(r) hW(r) 



(8) 



Having further decomposed each of the h[ fe l(r) matrices 
(k = 1 , . . . , 4) into NxN entries of 3 x 3 matrices as 
h|^(r) (I and m are counted along the row and the col- 
umn respectively), we have (with I as the identity matrix) 



•fi(r)Hl 



JT - (1 - e- a ° T ) 

a 



C > Sir: 



hfi(r) 



1 _ P — a r 
— 1 

ao 



+ J_ [ aor (i + e-aorj _ 2 + 2e - aoT } c \ 5 h , 



h^(T) = {7[l-e- Qor ]c} S lm and 



h£(r) 



-1-7 



t+ — (1 -e Qor ) 
a 



c \ Sim , (9) 



If we now form a 6A^ x 6A^ matrix G, [in the nota- 
tion of Eq. (§)] which looks like GR = = and 



-[2] 
3 lm 



-[31 
] lm 



where 



1 

1 
1 



(10) 



Then H(r) can be easily shown to satisfy p0|] 

[H(r)] T GH(r) = e- a " T G. (11) 

Thus, H(t) is generalized /i-symplectic with G, but it is 
not ^-symplectic, i.e., [H(r)] T J H(r) ^ e~ Q ° r J. The fact 
that the same analysis [Eqs. (|7|jTo|)] can be carried out for 
any constant coefficient of friction (not necessarily ao), 
implies that the CPR is exact for a collision-less gas of 
point particles obeying Eq. (^) with a constant coefficient 
of friction. This has been found previously in simulation 



data (TT|. 

For the transformation of over a binary collision 
between the i-th and the j-th sphere, we follow the ex- 
plicit derivation in Ref. Jl4| , which in turn is based on 
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the formalism developed simultaneously by Gaspard and 
Dorfman pl| , and by Dellago and co-workers p2| . The 
post-collisional infinitesimal deviation vector <5X + can be 
related to its pre-collisional value <5X_ by 



5X_, 



5X_ 



where the 6N x 6N matrix My can be decomposed into 
four 37V x 37V blocks, having the following structure 



= (I - 2NyNy) 



Here, R is a symmetric matrix. The above expression for 
M implies that the collisions are symplectic, i.e., 

M£- J My = J , 

but not generalized symplectic with G (M^GMy ^ G). 

We can now express the matrix L(t) in terms of the 
H and M matrices in the following way: if the dy- 
namics involves free flight segments separated by s in- 
stantaneous binary collisions at t±, t2,...,t s such that 
< ii < < 2 < • ■ ■ < t s < t, then 

L(i) = H(Ai s ) M lsJs H(At s _i) • • • M nn H(At ) . (12) 

Here, At s = t— t s and At; = tj+i— 1$ for i = 1, . .., (s— 1). 

The consequences of Eqs. (f^-p"2|) can be summarized 
by the following: for a collection of hard spheres obey- 
ing the SLLOD equations of motion with constant coef- 
ficient of friction oto, (a) the H matrices are generalized 
/i-symplectic with G, but not with J, and (b) the M ma- 
trices are symplectic but not generalized /i-symplectic 
with G. Hence, once the H and the M matrices arc com- 
bined together, as in Eq. (|l^), L(i) is seen to be general- 
ized /i-symplectic with neither G nor J. This is consistent 
with the claim that L(t) is not generalized /i-symplectic 
(and consequently, the CPR does not hold exactly) for a 
collection of hard spheres obeying the SLLOD equations 
of motion with constant coefficient of friction ao- 

The degree of deviation from an exact CPR must follow 
from the properties of L(i), and to estimate this devia- 
tion, we can use either K — G, or K = J in Eq. @. While 
the former choice implies that one has to try to estimate 
the deviation from an exact CPR from the distribution 
of the unit vectors N^'s and the collision angles for dif- 
ferent sets of binary collisions in the expression of M^-'s, 
the latter choice means that one can make the estimate 
by using the typical magnitude of a free flight time, i.e., 
the mean free time To. We choose the latter approach, 
because an estimate of the deviation from the exact CPR 
can be made at small 7, as a power series expansion in 
7. It is important to realize at this point that as the 
density sets a time scale in the form of the mean flight 
time To between collisions, the actual dimensionless small 
parameter corresponding to the shear rate is 7 = 7x0 . 



We begin by constructing another matrix Ho(Ai) by 
setting 7 = but ceo 7^ in the explicit form of H ( At) in 
Eqs. §|), i.e., 



Ho (At) = H(At)| 

QO^O, 7 — ■ 

It is easy to show that Ho (At) satisfies the equation 

[H (Ai)] T JH (Ai) = e~ aaAt J. 

Following Eq. ([l2|), we then form the matrix Lo(i) as 

L (t) = H (At,)M w , H (Ai s _!) • • • M hjl H {At ) , (13) 

such that all the My matrices in Eqs. ( |l2| ) and (13) are 
the same. Since both the My and the Ho(Ai) matri- 
ces are now /i-symplectic with J, so is \-o(t). As a con- 
sequence, the logarithms of the eigenvalues of \-o(t) = 
[Lo(£)] T Lo(£) pair exactly. This implies that if we arrange 
the corresponding Lyapunov spectrum 



Ao 



lim 

t— »oo 



Lo(«) 



l/(2t) 



in the decreasing order of magnitude as > X^' > 



>(0) 



> A, 



(0) 



l(0) , x(0) 



■v 6Ar , then A, ' + A 6 Xr_. i+1 = -ao- 
It is a simple exercise to show that H(A£) — Ho (At) = 
0(7 3 ), from which we conclude that for At = r = 0(tq) 

L(t) = L (t)[I +7 3 B], (14) 

where the matrix B is of order 1 in 7 and order 1 in N. 
Note that B contains higher powers of 7 as well. Be- 
cause it involves the matrix c and contributions from 
collisions between spheres, B is not proportional to I, 
and hence, we cannot regard it simply as a scalar fac- 
tor (in which case the exact conjugate pairing would be 
easy to obtain again). Equation ( |T^ ) implies that for 
AL(r) ^ L(r) - L (r) 



AL(r) 



7 3 [B T L (r) 



L (r)B] + 7 6 B T L (r)B. (15) 



From Eqs. ( |l4|) and (|l5|), we can now see that the dif- 
ferences between L(t) and Lo(r), and between L(r) and 
Lo(t) are small, by a relative order 7 3 . Therefore the 
logarithm of the eigenvalues of L(r) and Lo(r) also differ 
by a term of order 7 3 in an absolute sense. If we now di- 
vide the logarithms of these eigenvalues by the time r, we 
see that the finite time (for time r) Lyapunov exponents, 
calculated from Lo(r) and from L(r) (which we denote as 
A| (t) and \i(r) respectively, for i = 1, 2 . . . 6iV), differ 
by a term 0(7 3 /r) = 0(77 2 ). 

We make one further observation at this stage. The 
Lyapunov exponents (even the finite time ones) are in- 
variant under 7 — > —7, so in a power series expansion in 
7 p3], the odd powers vanish. Hence, we conclude that 
the logarithm of the eigenvalues of L(r) and Lo(r) must 
differ by a term of order 7 4 , i.e., the conjugate pairing of 
Ai(r)'s must be valid up to corrections of the form 77 s . 
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To explicitly extend this formalism to large t and 
thereby obtain a relation between A;S and A$ s, we need 
to sequentially concatenate a lot of L(r)'s. In general, 
these matrices neither commute with each other, nor with 
the B's, which prevents us from explicitly demonstrating 
how the deviation [L(t) — Lo(i)] is built up. However, we 
can argue in the following manner: L(i) and L (t) are pos- 
itive definite and symmetric. This allows us to express 
them in the form Lo(t) = exp(Ao) and L(t) — exp(A), 
where for large t, both the eigenvalues of Ao and A must 
behave ~ t. From this perspective, the difference be- 
tween the Lyapunov exponents for L(t) and Lo(£) is re- 
lated to (A — Aq). Since the difference between L(t) and 
\-o(t) has an explicit prefactor of 7 3 , so does A — A . 
Using the symmetry argument that the Lyapunov expo- 
nents have to be even functions of 7, we obtain 

A, + A 6JV - 4 +i = -a + 0(77 3 ) • i = l,...6JV (16) 

For the largest and the most negative Lyapunov expo- 
nents, it has been possible to show that they pair to — ao 
plus corrections of 0(77 3 ) by means of a kinetic theory 
approach ^4j,^5|, based on the independence of subse- 
quent collisions of a sphere. Likewise, one expects that 
in subsequent time-intervals of 0(tq), the L(r) matri- 
ces are not qualitatively much different from each other. 
Therefore, we expect that the coefficient of the 0(77 3 ) 
term in Eq. (pL|), to be of the same order as that for a 
flight time r = 0(tq) [i.e. of the order of B = 0(1)], and 
therefore Eq. ( |l6| ) to hold. 

In summary, for the SLLOD equations with a constant 
ao thermostat, the finite time Lyapunov exponents obey 
the CPR up to 0(77 3 ) when that time is of the order 
of the mean flight time, and this is expected to hold for 
the infinite time Lyapunov exponents too. Moreover, the 
iso-kinetic Gaussian thermostat is equivalent to the con- 
stant multiplier thermostat in the thermodynamic limit 
fl9f , and hence one expects that with an iso-kinetic Gaus- 
sian thermostat between collisions, the Lyapunov expo- 
nent spectrum also exhibits (9(77 3 ) deviations from the 
CPR, in the thermodynamic limit. Finally, given that 
the source of the CPR violation is basically the ao7CR 
term in Eq. (Q), one can argue that when the gas parti- 
cles interact with each other by means of a short-ranged, 
repulsive potential with a constant multiplier thermostat, 
the violation of the CPR would also be at least of 0(7 4 ) 
(for gas particles interacting with each other by means 
of a short-ranged, repulsive potential with an isokinetic 
Gaussian thermostat, the same results are expected) [ 15[ . 
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